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LONG-TERM  GOALS 

The  long-term  goal  of  this  project  is  to  develop  a  numerical  model  for  the  simulation  of  regional 
processes  characterized  by  a  disparity  of  temporal  and  spatial  scales,  such  as  nonlinear  internal  waves. 

OBJECTIVES 

To  develop  a  computational  model  of  large  amplitude  internal  gravity  waves  that  uses  adaptive  mesh 
refinement  in  combination  with  a  numerical  scheme  of  order  four  or  higher,  exhibiting  minimal 
numerical  dissipation  and  dispersion.  Furthermore,  to  have  the  capability  of  treating  realistic 
topography,  thus  requiring  the  capability  to  describe  complicated  geometry  and  also  parallel  execution 
to  reduce  the  large  computational  times  intrinsic  to  three-dimensional  simulations. 

APPROACH 

Physical  model 

We  write  the  equations  of  motion  as 

V  •  u  =  0 

p,  +(u  ■  V )p  -  N2(\  +  J3  p)v  =  0 
(\  +  ftp)\ut  +(u  -V )u\  =  -Vp  -(p  -  p)j  +—V2u+F 

p 

Here  all  quantities  are  assumed  to  nondimensionalized  with  reference  to  a  wave  length  scale  ^  ,  a  time 
scale  0  where  0  is  an  average  Brunt- Vaisala  frequency  and  an  average  background 
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maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1 .  REPORT  DATE  2.  REPORT  TYPE 

30  SEP  2007  Annual 

4.  TITLE  AND  SUBTITLE 

Holistic  Approach  To  Modeling  Non-Linear  Internal  Waves  With  A 
Three-  Dimensional  Nonhydrostatic  Model 

6.  AUTHOR(S) 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  North  Carolina,  Chapel  Hill, Chapel  Hill, NC, 27599 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

code  1  only 

14.  ABSTRACT 

15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF 

ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

5 

3.  DATES  COVERED 

00-00-2007  to  00-00-2007 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


2 

density.  We  introduce  the  Boussinesq  parameter  @  LN0I  g  wjt]1  g  gravitational  acceleration.  The 
reduced  density  P  is  introduced  such  that  @  P  P  is  the  density  perturbation  from  its  background 
stratification  given  by  P  Y  .  Similarly  a  reduced  pressure  is  introduced  such  that  P  P  is  the  pressure 
perturbation.  The  ^  term  captures  external  forces  such  as  Coriolis  forces  and 

=  P 
p0L2N0 

is  the  Reynolds  number.  In  the  Boussinesq  approximation  (  @  ^  )  we  set  ®  in  inertial  terms  but 

keep  a  finite  @  in  the  buoyancy  related  terms  (those  exhibiting  the  product  N  P  ). 

Numerical  treatment 

Domain  discretization.  The  domain  of  interest  is  simulated  on  logically  Cartesian  grids  in 
computational  space.  These  can  be  transformed  into  moderately  complicated  representations  of  real 
topography  by  employing  grid  mappings. 

Operator  splitting  approach.  The  equations  of  motion  are  advanced  in  time  by  successively 
considering  the  various  types  of  terms  that  appear.  The  hyperbolic  terms  resulting  from  advection 

P,  +(u-V)p  =  0 

ut  +  (u  •  V )u  =  —  V  « 

P 

are  updated  first  using  an  explicit  method  and  with  the  pressure  gradient  approximated  from  the 
previous  time  step.  Here  it  is  assumed  that  the  time  step  is  not  limited  by  instability  of  the  viscous 
tenn.  In  this  step,  the  quantities  that  are  advanced  in  time  are  defined  on  the  cell  centers,  whereas  the 
advective  field  is  defined  on  the  cell  edges. 

Then  the  effect  of  source  terms  (gravitational  acceleration  and  terms  arising  from  the  coordinate 
transfonnation)  are  added  to  the  velocities  defined  on  the  edges  of  the  computational  cells,  as  well  as  a 
suitable  interpolations  of  the  fluxes  calculated  for  the  cell  centered  velocities 

Ut=-(p-j3N2p)j+F. 

Finally  we  perform  a  projection  onto  the  space  of  divergence-free  velocity  fields  to  enforce  the 
constraint 

V  •  u  =  0 

on  the  edge  velocities.  The  potential  field  is  then  applied  to  the  cell  centered  velocities  as  well  to 
synchronize  the  two  discretized  versions  of  the  same  field.  This  operator  splitting  is  carried  out  at  each 
stage  of  a  low-dissipation  Runge-Kutta  time  stepping  procedure. 

Grid  adaptivity.  In  adaptive  mesh  refinement  one  must  identify  regions  requiring  additional 
resolution  by  evaluating  the  local  truncation  error.  Typically  this  is  done  through  a  trial  step.  This 
becomes  expensive  for  the  full  equations  of  motion  due  to  the  effort  involved  in  solving  the  Poisson 
equation  arising  from  enforcing  the  divergence-free  constraint.  To  avoid  this,  we  assume  that  the 
hyperbolic  terms  alone  give  a  sufficiently  accurate  indication  of  regions  requiring  additional  resolution 
and  let  the  grid  adaptivity  be  driven  by  trial  steps  of  this  stage  of  the  operator  splitting. 
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Poisson  and  Helmholtz  solvers. 


Each  subgrid  inherits  the  boundary  conditions  for  the  projection  scheme  from  the  edge  velocities  of  the 
parent  grid  already  projected,  so  that  mass  conservation  is  always  enforced.  To  setup  and  solve  the 
necessary  elliptic  solvers,  we  have  developed  a  package  that  breaks  down  the  problem  in  a  series  of 
steps:  The  linear  operators  corresponding  to  gradient  and  divergence  (in  curvilinear  coordinates)  are 
created  as  a  set  of  sparse  matrices  in  an  appropriate  representation;  next,  the  numerical  laplacian  is 
created  as  the  matrix  product  of  the  divergence  and  gradient  operators,  ensuring  consistency;  finally, 
the  resulting  linear  problem  Ax=b  is  solved  using  existing  highly  efficient  sparse  matrix  solvers.  The 
approach  is  object  oriented,  and  the  resulting  standalone  package  (ELLIKIT)  can  be  used  to  deal  with 
similar  problems  whenever  the  field  to  project  is  defined  on  cell  edges.  In  particular,  we  have 
implemented  hydrostatic  and  weakly  nonhydro  static  approximations  that  are  completely  transparent  to 
the  end  user. 

WORK  COMPLETED 

The  main  trust  this  year  was  devoted  to  implement  the  elliptic  and  parabolic  solvers  in  the  self- 
consistent  way  described  above,  in  generalized  coordinates.  We  also  considered  the  generation 
problem  of  internal  waves  along  the  NJ  shelf.  The  SW06  experiments  show  that  arrival  time  at  a  given 
location  across-shelf  of  NLIWs  is  not  well  predicted  by  the  baroclinic  tidal  phase  (unlike,  e.g.,  the 
South  China  Sea  case).  Furthermore  the  largest  NLIWs  appear  to  originate  during  neap-tide.  In  fact, 
the  barotropic  tidal  flow  is  too  weak  to  be  able  to  create  directly  NLIWs  by  the  standard  mechanism  of 
flow  over  topography.  In  order  to  investigate  other  possible  effects,  we  have  considered  what  effect  the 
presence  of  a  shelf-break  front  has  on  the  generation  process.  It  is  well  known  that  a  shelf-break  front 
is  present  in  the  area  year-around.  We  have  thus  conducted  a  series  of  experiments  in  which  different 
combinations  of  shelf-break  strength  and  stratification  are  considered. 

RESULTS 

Model:  Figure  1  shows  the  effect  of  the  shelf  break  front  on  the  development  of  the  internal  tide  as  the 
barotropic  tide  flows  over  the  shelf  edge.  Here  we  show  the  density  field  at  the  end  of  the  ebb  phase. 
On  the  left,  the  presence  of  the  front  causes  a  much  stronger  baroclinic  response  near  the  shelf  edge 
relative  to  the  control  case  without  front  (right).  In  both  cases,  the  barotropic  tide  is  subcritical. 

PUBLICATIONS 

A.  Scotti,  R.  Beardsley  and  B.  Butman,  “Generation  and  propagation  of  nonlinear  internal  waves  in 
Massachusetts  Bay”,  in  press  (2007). 

A.  Scotti,  J.  Wendelbo  and  S.  Mitran,  “A  note  on  the  parameterization  of  NLIWs  in  ocean  models”,  to 
be  submitted  to  Ocean  Modeling  (2007). 

Software 

The  implementation  of  the  above  method  are  to  be  made  publicly  available  to  the  wider  research 
community  in  accordance  with  the  grant  objectives.  Snapshots  are  periodically  made  available  in  the 

/  applications/2d/Boussinesq  and  /  applications/3d/Boussinesq  directories  of  the 
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BEARCLAW  package  (http://www.amath.unc.edu/Faculty/mitran/bearclaw.html).  ELLIKIT  is 
available  at  ftp://ftp.unc.edu/pub/marine/ascotti 

FURTHER  WORK 

On  the  modeling  side,  we  are  currently  exploring  how  to  parameterize  NLIWs  within  hydrostatic 
models  such  as  ROMS.  We  are  also  systematically  exploring  the  effects  of  the  shelf-break  front  on  the 
generation  of  NLIWs  during  SW06. 
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Figure.-l  [The  figure  shows  slices  of  the  density  field  near  the  shelf  edge.  One  plot 
shows  the  field  when  the  shelf-break  front  is  present,  the  other  without  a  shelfbreak 
front.  In  the  former  case,  the  baroclinic  response  is  much  stronger.] 
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